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Abstract 



Splitting methods for the numerical integration of differential equations 
of order greater than two involve necessarily negative coefficients. This 
order barrier can be overcome by considering complex coefficients with 
positive real part. In this work we review the composition technique used 
to construct methods of this class, propose new sixth-order integrators and 
analyze their main features on a pair of numerical examples, in particular 
how the errors are propagated along the evolution. 
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1 Introduction 

Splitting methods for the numerical integration of differential equations consti- 
tute an appropriate choice when the associated vector field can be decomposed 
into several pieces and each of them is explicitly integrable. 
Given the initial value problem 



can be integrated exactly, with solutions x(h) = ip^ (xo) at t = h, the time step. 
Splitting methods intend to approximate the exact flow iph by a composition of 
flows ip h . For instance, 
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both provide first-order approximations to the exact solution, since Xh(xo) = 
fh(xo) + 0(h?) and similarly for x*h (which is called the adjoint of Xh an d 
verifies x* h = XZ\)- 

It is possible to get higher order approximations by introducing more maps 
with additional real coefficients, <p a .. h , in ((3J). Perhaps the most popular split- 
ting method is the second order symmetric composition 

e [2] * [ml [2] [1] [2] [ml 

= Xh/2 O X h /2 = V h /2 ° • • • ° ° <Ph ° </2 ° • • • ° ( 4 ) 

When / in flT} is separable in two parts the above particularizes to 

[2] [1] * [1] [2] c [2] [2] [1] [2] 

Xh = <Ph °<Ph, Xh = <Ph°<Ph, S h =( Ph/2 0( Ph 0( Ph/2> ( 5 ) 

[21 

and S h is known as the Strang splitting |22] , the leapfrog or the Stdrmer-Verlet 
method [26], depending on the context where it is used. More generally, one 
may choose the coefficients a,, bi to achieve order r with the composition 

, [2] [1] [2] [2] [1] [2] fR , 

It turns out that iph can also be written in terms of Xh an d Xh 

I ( [2] [1] ^ o • • • o [ [2] o [1] ^ o ( [1] [2] 
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as long as 
Equivalently, 



1 Xa 2 fc Xaih (7) 
«2j-l + a2j, fy+i = a2j + «2j + l- (8) 



ai=h, a 2 j+\ = h + - a k ), a 2 j = ^(a fe - b k ), (9) 

k=l k=l 

with «o = «2s+i = 0. This relation remains valid if Ylt=i a i = Y2t=i ^ [IS]- ^ 
relevant consequence of this property is that, starting with the coefficients ai, bi 
of a given splitting method, we can get the coefficients ai for the composition 
([7]), which can be then applied in a more general setting with the maps Xh 
and x*h °f ©• A particular case widely used in practice to achieve high order 
approximations consists in considering compositions using the Strang splitting 
(dJ) as basic method, 

^ = s^ h o...os^ h os^ ih . (10) 

Splitting methods are, in general, explicit, easy to implement and pre- 
serve structural properties of the exact solution, thus conferring to the nu- 
merical scheme a qualitative superiority with respect to other standard inte- 
grators, especially when long time intervals are considered (see [6] for a re- 
view). Examples of these structural features are symplecticity, volume preser- 
vation, time-symmetry and conservation of first integrals. In this sense, split- 
ting methods constitute an important class of geometric numerical integrators 
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It has been shown that some of the coefficients in splitting schemes ([6]) are 
negative when the order r > 3 [91 [2T1 124] . In other words, the methods always 
involve stepping backwards in time. An elementary proof of this feature can 
be worked out as follows. It is quite straightforward to check that one of the 
necessary condition for the composition ([7]) (respectively, (|1U|) ) to have order 
r > 3 is 

k 
i=l 

with k = s (respect, k = 2s). Obviously, this sum vanishes only if at least one 
of the a.i is negative. In consequence, the flows tp^\ j = 1. . . . . m — 1 in (j4"|) and 
(pV 1 , j = 2, . . . , m — 1 in (|3J) evolve with at least one negative fractional time 
step. On the other hand, by taking into account the link (|8|) among coefficients 
of ([6]) and (J7J), condition (fTTj) with k = 2s leads to 

s 

^2( a 2j-i+ a 2j) = => ^ k / a 2k-i+ a lk < => a k = a 2 k-i+a 2 k < 0. 
i=i 

In a similar way, using the same condition with ao = a2s+i = 0, one has 

s 

^{alj + al j+1 ) = =► 3 I / a| + < ^ h = a 2 i + a 2 i+i < 

i=o 

and then at least one aj as well as one hi are negative. It must be stressed that 
condition (jlip still persists when the processing technique is used, so that the 
same conclusion also follows in this case [4]. 

In summary, the presence of negative coefficients in splitting methods of 
order higher than two is unavoidable if one restricts oneself to real coefficients. 
Of course, this does not suppose any special impediment when the flow of the 
ODE evolves in a group (such as in the Hamiltonian case), but may be unac- 
ceptable when the differential equation is defined in a semigroup [H], as occurs, 
for instance, with the simple heat equation ut = Au on the unit interval with 
homogeneous Dirichlet boundary conditions. Then the corresponding generated 
semigroup is well defined only for t > [llj. 

More generally, consider the nonlinear heat equation 

d d 

— u(x,t) = Y, D i( v i(x)D z u(x,t)) + F(x,u(x,t)) (12) 

i=l 

with functions V{ real and positive, and Di = d/dxi, on a certain domain 
n e R d . if a space discretization is carried out (either by finite differences or 
by a pseudospectral method), a large system of ODEs results which has to be 
numerically integrated in time. To this end, we can split the resulting equation 
into linear and nonlinear parts, but schemes of the form ([7]) or (|10j) of at most 
order r = 2 can only be applied, since the resulting discrete Laplacian with 
negative fractional time steps is not well conditioned. 
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(13) 



A technique used in practice to obtain the eigenvalues and eigenfunctions for 
a given potential V consists in numerically integrating the equation (after spa- 
tial discretization) along pure imaginary times (r = —it). Equivalently, the 
equation to be analyzed is 



which can be considered as a linear heat equation. The system evolves to the 
ground state whose norm decreases exponentially in proportion to the value 
of its energy (eigenvalue). By orthogonalization, one can make the system to 
evolve to any other eigenfunction [TJ [13] . In any case, whereas there is no special 
difficulties with numerically integrating equation (|13p using a splitting methods 
with negative fractional time steps, this is not the case for (|12p and (|14p due to 
the presence of the Laplacian. 

It has been noticed, however, that higher order splitting methods with com- 
plex coefficients having positive real part do exist [3j [16j [231 12H (25] ■ These 
schemes were reported mainly for theoretical purposes but received very lit- 
tle attention as practical numerical tools. Perhaps the main reason was that 
working with complex arithmetic makes the schemes more involved and, in 
many cases, also considerably more costly from a computational point of view 
(usually, four times more expensive). 

It is only within recent years that a systematic search for new methods 
with complex coefficients has been carried out and the resulting schemes have 
been tested in different settings: Hamiltonian systems in celestial mechanics 
[8], the time-dependent Schrodinger equation in quantum mechanics [19] 
and also in the more abstract setting of evolution equations with unbounded 
operators generating analytic semigroups [7JCEI]- In this sense, we recall that the 
propagator exp(zA) (z € C) associated with the Laplacian is well defined (in a 
reasonable distributional sense) if and only if Re(z) > [7J. More generally, it 
is possible to extend the semigroup related with parabolic PDEs into a sector 
in the right half plane of C [11] . 

The aim of this paper is to review some of the splitting methods with com- 
plex coefficients published in the literature, propose new sixth-order schemes in 
the class (110p and analyze them on a pair of simple numerical examples, to get 
a glance of the performance and main features of this kind of integrators and 
some of the difficulties involved. 

2 Integrators with complex coefficients 

Most of the existing splitting methods with complex coefficients have been 
constructed by applying the composition technique to the symmetric second- 
order leapfrog scheme S™. Thus, one gets a third-order method as 




(14) 




(15) 
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where the coefficients have to satisfy (jlip together with the consistency condi- 
tion 



a + /3 = 1 



2 ' 6 ' y 2 ^ 6 fi • 



Due to its simplicity, this scheme has been rediscovered several times, either as 
the composition (PT5j) [U [7] or by solving the order conditions required by 
© with s = 2 0QI]. 

A fourth-order integrator can be obtained with the symmetric composition 

cW _ o[2] «j[2] <j[2] nfil 

Although the necessary order conditions are the same, the time-symmetry of 
the composition rises the order by one (all the error terms at odd orders vanish 
identically) : 

2a + /3 = 1 \ _ _ 1 _ 2 x /3 e 2ifc7r /3 



2a A + f3 6 = J 2-2 1 /3 e 2'W3' ^ 2 -2 1 /3 e 2l W3 

with k = 0, 1, 2. Notice that for fc = 1, 2 it is true that Re(a), Re(/3) > 0. 

Another fourth-order method can be obtained by symmetrizing the third- 
order scheme (fT5|) . i.e., 

cW _ c[2] op] -[2] ~[2] , ^ 

Methods (fT5j) , (JT6J) and (fTT)) can be used to generate recursively higher order 
composition schemes as 

q [n+l] _ q [n] q [n] n& \ 

Here the coefficients have to verify the conditions a + /3 = 1, a n+1 + /5™ +1 = 0, 
whence 

_ 1 . sin(gj-Tr) J _«</<« _i if n is even, 

2 + '2 + 2cos(f±i7r) l-^^^ 1 if^isodd, 

and /3 = 1 — a. The choice I = gives the solutions with the smallest phase and 
allows one to build methods up to order six with coefficients having positive 
real part. This feature was stated in [25] and rediscovered in [TTj . 

In a similar way, one may use recursively a symmetric three term composi- 
tion, which allows to increase the order by two at each iteration: 

q [n+2] _ q [n] q [n] q [n] ,-. Q , 

with 

2a + (3 = 1, 2a n+1 + (3 n+1 = 0. 
The solutions providing coefficients with the smallest phase are 

e i7r/(n+l) 

a = 2 i/(n+i) _2e i7r /(™+ 1 )' / ? = 1 - 2 « 5 
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and methods up to order eight with coefficients having positive real part are 
possible. Moreover, methods up to order fourteen of the more general form (|10p 
with coefficients ctj with positive real part are attainable [3[TT]. An interesting 
(and open) question is to determine whether arbitrarily high orders can be 
attained or wether, as for the previous compositions, there is an order barrier 
for methods of the form (jlOp with Re(aj) > 0. Observe that any method of the 
form ()10p with coefficients having positive real part can be expressed in terms 
of the elementary flows tp^ with Re(A) > when S^ 2 ' is taken as the leapfrog 
([5]). This is also true, of course, in the more general case when / is split in m 
parts and S h is taken as the symmetric second order basic method ([!]) . 

For instance, suppose that / in (pQ) is separable in two parts, so that cSJ 2 ^ is 
given by ([5]). Then it is straightforward to check that the third order scheme 
(|15|) can be written as 

cP] I 2 ! W PI M I 2 ! / on \ 

S k = n; h Kzh </, nih ( 2 °) 

with a\ = \ + z^p, (ii = a\, b\ = ax/2, b% = 1/2, 63 = b\. This particular 
symmetry of the coefficients results in a method whose leading error terms at 
order 4 are all strictly imaginary [8]. 

Another question of practical nature is the construction of methods of the 
form (]10p with Re (ay) > involving the minimum number s of compositions 
for a prescribed order. For instance, the minimal number of compositions for 
achieving order 6 is s = 7. The corresponding order conditions can be written 
as [61 Q3H EH] 



2 °j = 

s 

where for each j = 1, . . . , s, 



a* = 0, fc = 3,5, (21) 
0, (k,£) G {(3, 1), (3, 2), (3, 3), (5, 1)}, (22) 



/ ,"7 



i-i 

a; 

a,; 



1=1 

This system of algebraic equations has several solutions with Re(ct,) > 0. 
Among them, we have chosen the two sets of coefficients collected in Table [TJ 
The first one corresponds to a symmetric method, a s+ ±-i = on, as scheme (fl~6j) . 
and was already found by Chambers [8]. The second method is apparently new, 
and possesses the special symmetry a s+ \-i = a*, as scheme (fT5|) (or (j20j) when 
expressed as p])). 



3 Numerical examples 

3.1 Example 1: the harmonic oscillator 

We consider the simple harmonic oscillator to illustrate some qualitative prop- 
erties of the previous composition methods with complex coefficients. That is, 
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Table 1: Coefficients of two 7-stage sixth-order methods of type (fTUj) : S76 is a 
symmetric method and S 7 6 is conjugate to a symmetric method (symmetric in 
the real part of the coefficients and skew-symmetric in the imaginary part). 

ai = 0.116900037554661284389 + 0.043428254616060341762 i 
a 2 = 0.12955910128208826275 - 0.12398961218809259330 i 
a 3 = 0.18653249281213381780 + 0.00310743071007267534 i 
a 4 = 0.13401673670223327014 + 0.15490785372391915239 i 

«5 = «3> "6 = "2, «7 = Oil 

S* 7 6 = 

at = 0.133741778914683628452 - 0.028839028371025553995 i 
a 2 = 0.12134019583938803504 + 0.11585180844272788007i 
a 3 = 0.13489797942731665044 - 0. 12906241362827633477 i 
a 4 = 0.22004009163722337213 



we take the Hamiltonian function H(q,p) = \{p 2 + q 2 ), with q,p € K. The 
corresponding equations of motion are linear and can be written as 

A B 

so that the numerical solution at time t = h furnished by method ([6]) is given 
by 

x{h) = K(h)x = e b °+ lhB e a ° hA e bshB ■ ■ ■ e b * hB e aihA e blhB x . (24) 

As is well known, for splitting methods with real coefficients the average error 
in energy remain constant for exponentially long times under suitable general 
conditions on the Hamiltonian. For the particular case of the harmonic oscil- 
lator and with a sufficiently small time step, this is true for all times, and the 
average error in positions grows only linearly. 

We propose here to check whether this also holds for methods with com- 
plex coefficients. To do that, we take as initial conditions (q,p) = (1,1) and 
integrate the system (|23|) for t £ [0, 200007r] using a constant time step. We 
measure the error in position and energy of the output obtained by propagating 
the solution with the splitting method and then computing the real parts of the 
results q ou t = Re(g), p ou t = Re(p). Figure [T] shows the results obtained with the 
following methods: (i) 523, the 2-stage third-order non-symmetric method fjl5|> . 
(ii) S34, the 3-stage fourth-order symmetric method (|16p . (iii) Sj6, the 7-stage 
sixth-order non-symmetric method, (iv) S76, the 7-stage sixth-order symmetric 
method. The coefficients of these two 6th-order methods are collected in Ta- 
ble [TJ The time step is chosen such that all methods require 27-28 evaluations 
per period. Notice the significant difference in the qualitative behavior of the 



7 



numerical solution. Whereas the error grows exponentially for the symmetric 
methods S3 4 and S76, this is not the case for S*23 and S76, which show a perfor- 
mance analogous to standard splitting methods with real coefficients: bounded 
energy error and linear growth of error in positions. Of course, such a behavior 
deserves a theoretical explanation, which we pursue next. 





Figure 1: Error in position and energy (taking the real part from the output) 
obtained with the 4th- and 6th-order symmetric schemes S34 and 576, and 
non-symmetric methods and S|6. The time step is chosen such that all 
schemes require 27-28 evaluations per period. 

The matrix K(h) in ()24|) is given explicitly by 

/ 1 W 1 a s h \ ( 1 ai h \ f 1 

m = { -b s+ ih 1 ) ( 1 )-( 1 ) ( -b lh 1 



In this way, one gets 

K(h) 



p(h)+d(h) q(h)+e(h) 
-q(h) + e(h) p(h) - d(h) 



where p(h), d(h) (respectively, q(h), e(h)) are even (resp. odd) polynomial 
functions having in general complex coefficients and det K(h) = p(h) 2 + q(h) 2 — 
d{h) 2 - e(h) 2 = 1. 

If the splitting method ([6]) is such that 



S 



a s - j+ i = a*, b s - j+2 = b* (25) 

(as happens, in particular, when it comes from a composition of the form (jlOp 
with a s -j + i = a*), then K(h)^ 1 = K(—h)*. More specifically, 

/ p(h) - d{h) -q(h) - e(h) \ ( p{hf + d{hf -q(h)* - e(h)* \ 
\q(h)-e(h) p(h)+d(h) J \q(h)*-e(h)* p(h)* - d(h)* J " 

This implies that p(y), q(h), and e(h) are real polynomials, whereas the coef- 
ficients of d(y) are purely imaginary. Notice that this is precisely the case of 
methods and S76. 

If, on the other hand, the splitting method is symmetric, i.e., it is of the 
form (J6]) satisfying 

ds-j+l = Oj; = bj 

(as happens, in particular, when it comes from a composition of the form (|10j) 
with a s —j+i = aj), then A'(/i) _1 = K(—h). This clearly implies that d(h) = 0, 
but in general the polynomials p(h), q(h), and e(h) have complex coefficients. 
For instance, methods <\16h (S^A) and 576 are such that p{h) has non-real coef- 
ficients. 

When a splitting method with matrix K{h) is used to integrate the harmonic 
oscillator, it is essential that p(h) € K. Otherwise K(h) n grows exponentially 
with the number n of steps. As a matter of fact, the eigenvalues of K(h) are 
Ai = e^W and A 2 = e'^, where 

<f>(h) = arccos(p(/i)), 

and thus max(|Ai|, | A 2 1 ) > 1 if p(h) £ M (and also if p(h) G R and \p(h)\ > 1). 
That is precisely the situation with methods S^A and S76, and thus they are 
useless when integrating harmonic oscillators or systems that can be considered 
as close perturbations of harmonic oscillators with the partition (j23|) . 

From the previous comments, it is clear that instability will take place when 
integrating the harmonic oscillator unless —1 < p(h) < 1. In fact, the numerical 
solution can still be (weakly) unstable when p(h) 2 = 1 with q(h)d(h)e(h) ^ 
[5]. Furthermore, it is shown in [5] that, for stable numerical solutions (that is, 
either —1 < p(h) < 1 or p(h) 2 = 1 with q(h) = d(h) = e(h) = 0), one has 

\ —sm(n<p(n)) cos(n0(n)) J 

with a suitable 2x2 matrix Q(h) (typically close to the identity matrix). In 
consequence, the numerical solution x n = (q n ,Pn) is such that x n := Q(h)x n 
corresponds to the exact solution at t n = nh of a harmonic oscillator with 
frequency u> = l/h(p(h) ~ 1. This feature explains why schemes and 
£76, when applied to the harmonic oscillator f)23|) with h = tt/7 and h = tv/2 
respectively, exhibit a linear error growth in positions and a bounded error 
in energy, since for such methods, p(h) = 1 — h 2 /2 + • • • is real and satisfies 
p{h) G (—1, 1) for the values of h considered in the numerical experiments. 
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3.2 Example 2: The Volterra-Lotka problem 

Consider now the Volterra-Lotka problem 

u = u(v-2), v = v{l-u). (26) 

This is a very simple nonlinear system which allows us to make a preliminary 
study about the behavior and performance of methods with complex coefficients 
in the transition process from a linear to a nonlinear problem. In a neighborhood 
of the steady state at (u*,v*) = (1,2) the system can be considered close to a 
harmonic oscillator. The nonlinear contributions are manifest as we move away 
from it. The system evolves along periodic trajectories around the equilibrium 
point in the region < u, v determined by the first integral I(u, v) = ln(uv 2 ) — 
(u + v). 

The vector field f(u, v) = (u(v—2),v(l—u)) can be separated in two solvable 
parts and this can be done in different ways. We consider the following split: 
/a = (u(v — 2),0) and fs = (0,^(1 — u)) (although the linear and nonlinear 
separation can also be considered). 

We take as initial conditions (uq, vq) = (2, 4), integrate up to t = 20000 x 2ir 
and measure the relative error in the first integral, \I — Io\/\Io\. As in the 
previous example, we integrate using complex arithmetic and take the real 
parts of u and v only for representing the output. Figures [2]-( a) and (b) show 
the results obtained for time steps h = 4^ and four times smaller h = 
with m the number of stages of each method. In this way, all methods require 
the same number of evaluations. Contrarily to the pure harmonic oscillator, 
we observe a secular error growth in the determination of the first integral 
for all methods which diminishes considerably when the time step is reduced. 
The observed behavior resembles what takes place with the so-called pseudo- 
symplectic methods (integrators of order n which preserve symplecticity up to 
order p > n), where the dominant errors behave as £/ = Ch n + tDh p for some 
constants C and D. If p > n the secular part of the error does not manifest for 
relatively long times when the time step is reduced. 

We have repeated the same experiment but, after each time step, we discard 
the imaginary part of u and v and initiate the next step only with their real 
part. In other words, we project each component on the real axis at the end of 
each integration step. The results obtained are shown in Figures [2]-(c) and (d). 
Obviously, this way of proceeding does not preserve symplecticity any more but 
the results obtained suggest that a significant improvement in accuracy can be 
achieved. 

4 Conclusions and outlook 

We have presented a short review of the splitting and composition technique to 
build methods of order greater than two with complex coefficients with positive 
real part. This procedure allows to overcome the order barrier where splitting 
methods of order greater than two involve necessarily negative coefficients in 
the real space. In general, splitting methods with complex coefficients are con- 
siderably more expensive than the corresponding methods with real coefficients 
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Figure 2: Relative error in the first integral / = \n{uv 2 ) 
Volterra-Lotka problem with initial conditions (uq,vq) = 



steps h 



fff and h 



— (u + v) for the 
(2, 4) for the time 



210' 



with m the number of stages of each method. 



(about four times more expensive), and this make them hardly competitive in 
practice. For this reason, one can think that the main application of the new 
methods could be on parabolic PDEs, where higher order methods with real co- 
efficients (which necessarily have some negative coefficientes) can not be used. 
However, there is a number of problems which evolve in the complex space 
where using methods which complex coefficients does not necessarily mean in- 
creasing the cost of the algorithm. This can be the case, for instance, of the 
Schodinger equation (fl~3|) . 

As for the practical implementation of splitting methods with complex co- 
efficients, in [8] it is claimed that one has to carry the numerical integration in 
complex variables, and (for problems with real solutions) one should take either 
the real part of the variables or their modulus only for the output. However, we 
have observed that removing the imaginary part at each step, i.e. projecting on 
the real space at each step, the error grow can be considerably diminished in 
some cases. In the numerical examples considered in previous section, the linear 
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error grow in the first integrals originate from different sources depending on 
wether the projection onto the real domain is performed after each step or not. 
In the first case, the projection after each step destroys symplecticity but only 
at a higher order, and the schemes can be considered as pseudosymplectic. In 
the second case, the method is actually symplectic and can thus be (formally) 
considered as an exact solution of a Hamiltonian system in the complex domain, 
which have qualitatively different properties to trajectories in the real domain. 
We have also noticed that the higher order methods present a considerably re- 
duced error grow. Then, it seems appropriate to look for efficient higher order 
methods with complex coefficients. In general, symmetric splitting methods are 
desirable. However, we have shown that for the harmonic oscillator symmetric 
methods (with non-real stability polynomial) present an exponential error grow, 
which is not the case for methods with the special symmetry (|25|) . In a prelimi- 
nary search of methods, we have presented a new sixth-order method with that 
special symmetry. This is an interesting subject to be further explored since 
many problems in different applications can be considered as perturbations to 
the harmonic oscillator. 
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